#include "Dwarf.Unmanaged.Omp.h"


Solver :: Solver()
{

}

Solver :: ~Solver()
{    
    delete hmm;
    delete vit;
}

//For the demonstration purposes of the project,
//we have chosen Viterbi algorithm for Hidden Markov Model
void Solver::solve() 
{

    /* 1. Initialization  delta and psi */
    int numberOfStates = hmm->getNumberOfStates();
    #pragma omp parallel for
    for (int i = 0; i < numberOfStates; i++) 
    {
        vit->delta[0][i] = hmm->initialStates[i] * (hmm->observationProbabilityMatrix[i][vit->observationSequence[0]]);
        vit->psi[0][i] = 0;
    }

    /* 2. Compute */
    int lengthOfObsSeq = vit->getLengthOfObservationSequence();
    for (int t = 1; t < lengthOfObsSeq; t++) 
    {
        //parallelization of inner cycle
        #pragma omp parallel for
        for (int j = 0; j < numberOfStates; j++) 
        {
            double maxval = 0.0;               //temporary storage for max value for delta array
            int maxvalind = 0;                 //temporary storage for max value for psi array
            double val = 0.0;
            //finding the maximum
            for (int i = 0; i < numberOfStates; i++) 
            {
                val = vit->delta[t-1][i]*(hmm->stateTransitionMatrix[i][j]);
                if (val > maxval) 
                {
                    maxval = val;
                    maxvalind = i;
                }
            }
            vit->delta[t][j] = maxval*(hmm->observationProbabilityMatrix[j][vit->observationSequence[t]]);
            vit->psi[t][j] = maxvalind; 
        }
    } 

    /* 3. Termination */

    vit->setProbability(0.0);

    vit->sequenceOfState[lengthOfObsSeq - 1] = 0;
    #pragma omp parallel for
    for (int i = 0; i < numberOfStates; i++) 
    {
        if (vit->delta[lengthOfObsSeq - 1][i] > vit->getProbability()) 
        {
            vit->setProbability(vit->delta[lengthOfObsSeq - 1][i]);
            vit->sequenceOfState[lengthOfObsSeq - 1] = i;
        }
    } 

    /* 4. Path (state sequence) backtracking */
    for (int t = lengthOfObsSeq - 2; t >= 0; t--)
    {
        vit->sequenceOfState[t] = vit->psi[t+1][vit->sequenceOfState[t+1]];
    }
}
